DSWA-TR-98-62 

AN  INNOVATIVE  COMPUTATIONAL  METHOD 
FOR  PREDICTING  NEUTRON  TRANSPORT 


Approved  for  public  release;  distribution  is  unlimited. 


Prepared  for: 

Defense  Threat  Reduction  Agency 
45045  Aviation  Drive 
Dulles,  VA  20166-7517 

DNA001-96-C-0103 


Jay  I.  Frankel 


Prepared  by:  University  of  Tennessee  -  Knoxville 

College  of  Engineering 
414  Dougherty  Engineering  Building 


20000828  006 


Knoxville,  TN  37996-2210 


J)TIO  £*?!/:  U 


DESTRUCTION  NOTICE: 


Destroy  this  report  when  it  is  no  longer  needed.  Do 
not  return  to  sender. 


PLEASE  NOTIFY  THE  DEFENSE  THREAT 
REDUCTION  AGENCY,  ATTN:  ADM,  6801 
TELEGRAPH  ROAD,  ALEXANDRIA,  VA.  IF 
YOUR  ADDRESS  IS  INCORRECT,  IF  YOU  WISH 
IT  DELETED  FROM  THE  DISTRIBUTION  LIST, 
OR  IF  THE  ADDRESSEE  IS  NO  LONGER 
EMPLOYED  BY  YOUR  ORGANIZATION. 


CUT  HERE  AND  RETURN 


I 


DISTRIBUTION  LIST  UPDATE 

This  mailer  is  provided  to  enable  DTRA  to  maintain  current  distribution  lists  for  reports.  (We  would 
appreciate  vou  providing  the  requested  information.) 

□  Add  the  individual  listed  to  your  distribution  list. 

□  Delete  the  cited  organization/individual. 

□  Change  of  address. 

NAME: _ 

ORGANIZATION: _ 

OLD  ADDRESS  NEW  ADDRESS 


Note: 

Please  return  the  mailing  label  from  the 
document  so  that  any  additions,  changes, 
corrections  or  deletions  can  be  made  easily. 
For  distribution  cancellation  or  more 
information  call  DTRA/ADM  (703)  325-1036. 


TELEPHONE  NUMBER:  (  ) _ 

DTRA  PUBLICATION  NUMBER/TITLE  CHANGES/DELETIONS/ADDITIONS,  etc.) 

(Attach  Sheet  if  more  Space  is  Required) 


DTRA  or  other  GOVERNMENT  CONTRACT  NUMBER: _ 

CERTIFICATION  of  NEED-TO-KNOW  BY  GOVERNMENT  SPONSOR  (if  other  than  DTRA): 
SPONSORING  ORGANIZATION: 

CONTRACTING  OFFICER  or  REPRESENTATIVE: 


SIGNATURE: 


DEFENSE  THREAT  REDUCTION  AGENCY 
ATTN:  ADM 

45045  AVIATION  DRIVE 
DULLES,  VA  20166-7517 


DEFENSE  THREAT  REDUCTION  AGENCY 
ATTN:  ADM 

6801  TELEGRAPH  ROAD 
ALEXANDRIA,  VA  22310-3398 


REPORT  DOCUMENTATION  PAGE 


form  Approved 
OMB  No .  0704-0188 


Public  reporting  burden  for  :ms  a>fl*ction  ot  .nfo/mjtton  .3  osiwnatod  to  ivto^  1  hour  p«r  roooo^s*.  including  mo  Tim*  ‘ot  r owowmg  .nstruci»ons.  soarr-rx;  cutting  d*t*  sourcos. 
gather  mg  and  maintaining  tho  data  noodod.  and  completing  and  reviewing  sr-  collection  of  information  Send  commern  rygordmg  this  burden  estimate  or  jny  otner  asoect  of  thn 
collection  of  information,  including  suggestions  for  reducing  this  burden,  to  Washington  Headouarters  Services.  Oirecmnrt  for  information  0  pent  tons  and  ;«oom.  !2iS  Jefferson 
Oaws  Highway.  Sunt  120*.  Arlmgton.  VA  22202-4302.  and  to  the  Office  o >  Management  and  fludget.  Piperworit  Reduction  a-nect  <0704-0188).  Wasneigton.  DC  :0S03 _ 


2.  REPORT  DATE  3.  REPORT  TYPE  AND  DATES  COVERED 

March'  2000  ! 


1 .  AGENCY  USE  ONLY  (Leave  blank) 


4.  TITLE  AND  SUBTITLE 

An  Innovative  Computational  Method  for  Predicting 
Neutron  Transport 


6.  AUTHOR(S) 

Jay  I.  Frankel 


Technical 


5.  FUNDING  NUMSERS 

C  -DNA  001-96-C-0103 
PE  -62715H 
PR  -AP 
TA  -D 

WU  -DH53052 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADORcSS(ES)  8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

University  of  Tennessee-Knoxville 
College  of  Engineering 
414  Dougherty  Engineering  Building 
Knoxville,  Tennessee  37996-2210 


3.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES)  10.  SPONSORiN 

Defense  Threat  Reduction  Agency  agency r_ 

45045  Aviation  Drive  DSWA-TR- 

Dulles ,  VA  20166-7517 

CPWP/Kehlet  _ 


1 1 .  SUPPLEMENTARY  NOTES 

This  work  was  sponsored  by  the  Defense  Threat  Reduction  Agency  under 
RDT&E  RMC  code  B  4662  D  AD  53052  3200  A  AF  25904D. 


SPONSORING  MONITORING 
AGENCY  REPORT  NUMBER 

DSWA-TR-98-62 


12a.  DISTRIBUTION/ AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  is  unlimited. 


12b.  DISTRIBUTION  C00E 


1 3.  ABSTRACT  (Maximum  200  words) 


This  final  report  indicates  the  work  performed  under  DNA-00 1-96-0 103  during  the  span  covering  June  1996 
through  June  1998.  Several  important  findings  are  reported  strongly  suggesting  that  the  proposed  numerical 
methodology  can  perform  as  initially  conceived  upon  resolving  the  remaining  technical  issue  involving  the 
numerical  integration  of  integral  containing  Hadamard  type  kernels.  Preliminary  calculations  are  presented  in 
both  cartesian  and  spherical  geometries  indicating  that  the  symbolically  augmented  method  has  merit  when 
addressing  the  numerical  simulation  of  the  integral  form  of  the  linearized  Boltzmann  transport  equation. 


U.  SUBJECT  TERMS 


Symbolic  manipulation,  Fredholm  integral  equations,  singular 
kernels.  Linearized  transport  equation 


15.  NUMBER  0?  RAGES 

56 


15.  PRICE  CODE 


17  SECURITY  CLASSIFICATION 
OF  REPORT 

18.  SECURITY  CLASSIFICATION 
OF  THIS  PAGE 

19.  SECURITY  CLASSIFICATION 
OF  ABSTRACT 

20.  LIMITATION  CF 

ABSTRACT 

UNCLASSIFIED 

UNCLASSIFIED 

UNCLASSIFIED 

SAR 

NSN  7S40-0l-2aO.Sfifin 

Standard  2E3  .Rev  2-39) 

Table  of  Contents 


Section  Page 

Figures .  iii 

Tables .  iv 

1  Summary .  1 

2  Concept  Development .  2 

2.1  Introduction .  2 

2.2  Symbolic  Form  of  the  Kernels .  2 

2.3  Orthogonal  Collocation  with  Chebyshev  Basis  (Isotropic  Scattering  Case) .  4 

2.4  Orthogonal  Collocation  with  Chebyshev  Basis  (Linear  Anisotropic  Case) .  7 

2.5  Orthogonal  Collocation  with  Chebyshev  Basis  (Anisotropic  Case  [12]) .  10 

3  Symbolic  Manipulation  of  the  Kernel  Functions .  19 

3.1  Identification  of  the  Kernel  in  Terms  of  Known  Functions .  19 

3.2  General  Anisotropic  Scattering  Cases .  19 

3.3  Identification  of  A'm,n(7 hO  m  Terms  of  Known  Functions .  20 

3.4 "Modified  Kernel  Approach . . .  22 

3.5  Identified  Anisotropic  Scattering  Test  Phase  Function .  23 

4  Conclusions  and  Recommendations .  25 

5  References .  26 

Appendix 

DAHC35-95-P-3171  [4]  .  A-l 


ii 


Figures 


Figure  Page 

1  $0(77)  and  $1(77)  distributions  for  various  values 

of  N  when  u>  —  0.5,  ai  =  —2.7  and  R  =  1 .  17 

2  $o(7?)  and  $1(77)  distributions  for  various  values 

of  N  whenu;  =  0.5,  ai  =  —2.7  and  R  =  5 .  18 

3  Simple  back-scattering  phase  function . . .  24 

4  Highly  lobular  phase  function .  24 


Tables 


Table  Page 

1  Numerical  results  for  the  reflectivity  using  the  proposed  Chebyshev 

collocation  method  for  various  u  and  R .  7 

2  Comparison  of  results  for  the  emissivity,  e  using  the  proposed  method 

with  that  of  a  Galerkin  solution  [12]  when  R  =  1  for  various  ui  and  N .  13 

3  Comparison  of  results  for  the  emissivity,  e  using  the  proposed  method 

with  that  of  a  Galerkin  solution  [12]  when  R  =  5  for  various  u>  and  N .  14 

4  Expansion  coefficients  for  $o  (v)  when  R  =  1,  oi  =  —2.7,  and  u>  =  0.-5 .  14 

5  Expansion  coefficients  for  $1(77)  when  R  —  1,  ai  =  —2.7,  and  uj  =  0.5 . 15 

6  Expansion  coefficients  for  $0  {v)  when  R  =  5,  ai  =  —2.7,  and  u  =  0.5 . 15 

7  Expansion  coefficients  for  $1(77)  when  R  —  5,  ai  =  —2.7,  and  u  —  0.5 . 16 


IV 


Section  1 
Summary 


This  final  report  indicates  the  work  performed  under  DNA-001-96-0103  during  the  span 
covering  June  1996  through  June  1998.  Several  important  findings  are  reported  strongly 
suggesting  that  the  proposed  numerical  methodology  can  perform  as  initially  conceived 
upon  resolving  the  remaining  technical  issue  involving  the  numerical  integration  of  inte¬ 
gral  containing  Hadamard  type  kernels.  Preliminary  calculations  are  presented  in  both 
cartesian  and  spherical  geometries  indicating  that  the  symbolically  augmented  method 
has  merit  when  addressing  the  numerical  simulation  of  the  integral  form  of  the  linearized 
Boltzmann  transport  equation. 


1 


2.1  Introduction. 


Section  2 

Concept  Development 


In  a  series  of  articles  [1-3],  a  symbolically  enhanced  methodology  has  been  developed  for 
solving  the  integral  form  of  the  Boltzmann  transport  equation  in  the  presence  of  anisotropy. 
The  integral  for  the  transport  equation  has  been  utilized  in  these  past  one-dimensional 
studies  since  a  reduction  in  dimensionality  occurs.  The  resulting  system  of  weakly  singular 
Fredholm  integral  equations  of  the  second  kind  involves  the  various  Legendre  moments  of 
intensity  (radiative  transport)  or  neutron  flux  (neutron  transport). 

In  [1],  an  important  computational  attribute  was  reported.  It  was  observed  that  symbolic 
could  be  used  to  automate  the  identification  of  the  kernel  function  needed  to  arrive  at 
the  integral  form  of  the  transport  equation  in  a  slab  geometry.  Theoretical  considerations 
involving  convergence  and  error  analyses  were  established  [2,3]  based  on  a  weighted  resid¬ 
ual  methodology  where  Chebyshev  polynomials  were  used  as  the  basis. set  in  the  series 
expansion  of  the  unknown  dependent  variables.  The  convergence  rate  was  established  for 
the  isotropic  scatering  case  when  implementing  orthogonal  collocation.  In  [3],  LaClair  and 
Frankel  investigated  linear  anisotropic  scattering  and  developed  rigorous  error  estimates 
based  on  functional  analysis  that  extended  the  error  estimates  developed  in  [2]. 

This  report  focuses  on  the  spherical  coordinate  problem  and  extends  the  numerical  findings 
reported  in  [4]  which  involved  highly  anisotropic  scattering  in  a  plane-parallel  geometry. 
This  singular  geometry  produces  additional  analytical  and  numerical  difficulties  that  do 
not  appear  in  the  plane-parallel  geometry  associated  with  a  one- dimensional  slab  geometry. 
Attention  is  first  directed  toward  the  identification  of  a  method  for  symbolically  represent¬ 
ing  the  kernel  function  in  the  spherical  coordinate  system.  Finally,  some  discussion  is 
directed  toward  developing  a  numerical  solution  to  a  series  of  coupled  integral  equations 
containing  these  newly  formed  kernel  functions. 

The  first  two-thirds  of  the  research  effort  was  highly  successful.  The  final  one-third,  which 
involved  numerical  integration  of  singular  integrand,  was  not  as  successful.  If  the  program 
would  to  continue  for  an  additional  year,  this  final  hurdle  should  be  resolvable. 

2.2  Symbolic  Form  for  the  Kernels. 

Unlike  Phase  I  of  this  investigation,  which  involved  the  integral  form  linearized  Boltzmann 
transport  equation  in  a  plane-parallel  geometry,  symbolic  implementation  in  the  spheri¬ 
cal  geometry  is  not  as  straightforward.  This  was  seen  during  July  1996.  Fortunately,  an 
observation  was  made  making  use  of  a  new  intermediate  form  which  allows  for  the  sym¬ 
bolic  implementation  to  be  accomplished.  This  observation  does  not  follow  conventional 
expressions  used  in  describing  the  kernel  functions. 

As  noted  by  Thynell  and  Ozisik  [5],  the  necessary  kernel  for  the  integral  form  of  the 
transport  equation  in  the  spherical  coordinate  system  is 
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where  Pm(u)  represents  the  mth  Legendre  polynomial  of  the  first  kind.  Letting  y  =  r  +  x 
and  z  =  \r  —  x\  allows  Eq.  (l)<to,be  given  as 
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Unlike  the  slab  geometry  [4] ,  some  additional  algebraic  effort  and  some  new  function  defini¬ 
tions  must  be  introduced  in  order  to  arrive  at  a  form  conducive  to  symbolic  manipulation. 
For  illustrative  purposes,  we  shall  consider  the  case  where  m  —  n  =  1,  therefore  our 
attention  is  focused  on 


Kn(r,x) 


-x2  +  t2 
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N, 
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where  Pi  ( u )  =  u.  Expanding  the  polynomial  product  of  the  Legendre  polynomials  of  the 
first  kind  and  then  simplifying  and  collecting  terms  produces 
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In  general,  we  can  define 
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and  use  this  definition  in  developing  a  general  expression  for  the  kernal  function.  Using 
Eq.  (7),  the  integral  can  be  analytically  evaluated  to 

Gj(y,z)  =  -^Ej(z)--^E1(y),  j=  0,1 .  (8) 

This  definition  for  Gj(y,z )  can  be  used  to  describe  any  degree  of  anisotropy.  For  example, 

Koo(r,x)  =  Gi(y,z),  (9) 

K0i{r,x)  =  -  2r---- G2{y,z )  +  ^  J^_  €_tdt,  (10) 

Kio(r,x)  =  2x~G^V'  Z^~  2x  Jt-  (n) 
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“  in  °3iy '  2>  -  4^ /_2  (12) 

where  y  =  r  +  x  and  z  =  |r  -  ar|.  Also,  the  remaining  integrals  displayed  in  Eqs.  (10-12) 
can  clearly  be  integrated  analytically. 

This  procedure  can  continue,  though  as  m  and  n  increase,  the  algebraic  level  of  difficulty 
also  incfeases.  The  functions  described  in  Eqs.  (9-12)  have  been  verified  using  the  study 
performed  by  Abulwafa  [6]. 


2.3  Orthogonal  Collocation  with  Chebyshev  Basis  (Isotropic  Scattering  Case). 

A  preliminary  simulation  was  performed  to  indicate  the  accuracy,  relative  speed,  and 
memory  requirements  for  the  one-speed,  bare  sphere,  isotropic  scattering  albedo  problem 
[2,5,6-10].  This  problem  is  well-studied  and  the  open  literature  offers  benchmark  numerical 
results  as  obtained  by  the  F-N  method  [8]  and  by  the  Galerkin  method  with  a  nonorthog- 
onal  basis  set  [7].  These  results  report  the  numerical  value  for  the  albedo  for  various 
scattering  coefficients  and  radii.  Also,  both  [8]  and  [7]  report  the  number  of  terms  neces¬ 
sary  in  the  series  representation  to  obtain  five  digit  accuracy  for  the  albedo.  It  should  be 
noted  that  the  calculation  of  the  albedo  requires  fewer  terms  to  achieve  five  digit  accuracy 
than  the  actual  neutron  flux.  This  occurs  because  of  the  integral  definition  of  the  albedo 
in  terms  of  the  zeroth  Legendre  moment  of  the  neutron  flux. 

In  order  to  verify  and  evaluate  the  accuracy  of  the  approach,  preliminary  calculations  are 
developed  using  the  radiative  heat  transfer  community  results.  Tabular  results  have  been 
reported  by  Ozisik  and  his  co-workers  over  the  last  decade. 

Defining  the  mth  Legendre  moment  of  the  intensity,  I (r,  p)  as 

Gm(r)  =  f  Pm(p)I(r,p,)dp,  m  =  0, 1, ...,  N,  (13) 

Jn=- 1 

where  Pm(z)  represents  the  mth  Legendre  polynomial  of  the  first  kind.  With  this  defini¬ 
tion,  the  linearized  Boltzmann  transport  equation  can  be  expressed  in  an  integral  form. 
Following  Thynell  and  Ozisik  [5],  the  general  integral  form  of  the  transport  equation  in 
spherical  coordinates  can  be  expressed  as 


R  N 

rGn(r)  =  2tt  J _  \^Q{x)KQn{r,x)  +  ~  bmGm{x)Kmn{r,x ) 

x  ^  m=0 


xdx+ 
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where  the  kernel  function  Kmn(r,x)  is  given  by 
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Kmn(r,  X)  =  Pm  (  —  )  Pn  (  —  )  —dt,  771,  71  =  0,  1, 
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while 

o(x, r,  p)  =  \J  x2  —  r2(l  -  /x2), 

(16) 

and 

Po(r,ft)  =  \/l  ~  {r/R)2(l  -  fj,2), 

(17) 

Assuming  a  transparent  boundary  at  r  =  R  allows  us  to  write  the  boundary  condition  [7] 

I~(R,  -p)  =  f(p),  (18) 

In  deriving  Eq.  (17),  the  conventional  form  of  the  phase  function  was  assumed;  namely, 

N 

p(/X,/i')  =  bmPMPm(p'),  bo  =  1,  (19) 

771= 0 

where  the  {fem}^=0  are  known  coefficients. 

In  the  case  of  isotropic  scatter  (IV  =  0),  the  kernel  displayed  in  Eq.  (15)  reduces  to 

Koo{r,x)  =  Ei{\r-x\)-Ei(r  +  x),  z,re[0,.R],  (20) 

where  E\{z)  represents  the  first  exponential  integral  function  which  contains  a  known 
logarithmic  singularity  as  z  — *  0. 

In  order  to  introduce  the  Chebyshev  basis  functions  later,  we  need  map  the  domain  from 
re[0, 1?]  to  77e[ — 1, 1] .  Clearly, 

r  =  A(1  +rj),  re[0  ,R],  r,e[-l,  1],  _  (21) 

where  A  =  Rj 2  accomplishes  the  required  transformation.  Using  Eq.  (21),  the  transport 
equation  shown  in  Eq.  (14)  can  be  written  as 

f,\2  r 1 

A(1  +  77) (77)  =  A(1  +  v)P(v)  ■+ — 2~  J  (1  +  OViOkooiViOdi,  77c  [-1,1],  (22) 

where 

^(77)  =  Go(A(l  +  77)),  (23) 

kooiv,  0  =  AToo(A(l  +  77),  A(1  +  £))  =  Ei(X\rj  -  Cl)  -  Ei{\{2  +  rj  +  0),  (24) 

while  ^ 

F(77)  =  4tt/  cosh[A(l  +  7 ?)p]e- Ve[-l,l],  (25) 

Jn= 0 

where  f{p)  —  1,  i.e.,  isotropically  incident  source  at  r  =  R. 

Let  the  solution  for  $(77)  be  expressible  in  terms  of  an  infinite  series;  namely, 

OO 

^(v)  =  ^2ajTj(/n),  *?<-!,  1],  (26) 

3— 0 
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where  {Tj(r})}jL0  represent  Chebyshev  polynomials  of  the  first  kind.  The  unknown  ex¬ 
pansion  coefficients  are  denoted  by  the  infinite  set  of  constants  practice,  this 

infinite  series  representation  must  be  truncated  after  a  finite  number  of  terms,  say,  N  +  1. 
Therefore,  we  form  an  approximation  to  (77)  which  is  denoted  as  $^(77)  an<^  is  given  by 

N 

'  #(77)  «  ^(77)  =  77e[-l,l].  (27) 

3— 0 

Upon  subtstituting  the  assume  expansion  for  the  unknown  function  ^ s{r])  into  Eq.  (22), 
the  residual  equation  is  formed;  namely, 

N  U)\2 

Rn{^n{v))  =  A(1  +v)F(v)  ~  [A(l  +  77)^/(77)  ~  -5 ^[-1,1],  (28) 

3= 0 

where 


Aj(r])  =  [  (l  +  OTj(Okoo{v,Z)<%,  j  =  0,1,...,  N,  T7e[— 1,1],  (29) 

Ji=-i 

The  function  defined  in  Eq.  (29)  can  be  integrated  analytically  without  difficulty.  Numeri¬ 
cal  evaluation  of  the  weakly  singular  integrand  can  also  be  accomplished  by  first  performing 
one  or  more  integration  by  parts  prior  to  applying  a  numerical  quadrature. 

For  illustration  and  comparison  purposes,  we  propose  to  use  orthogonal  collocation.  Un¬ 
like  the  Galerkin  method,  this  method  does  not  require  any  additional  integrations  to 
be  performed.  The  collocation  method  is  defined  through  the  following  inner  product 
(orthogonality)  statement 

^RN(^N(v)),S(v-Vk)^  =0,  k  =  0,l,...,N.  (30) 

Here,  the  collocation  points  are  denoted  as  77*.  k  =  0, 1, ...,  N  and  given  by  the  open  rule 


T]k  —  cos 


f(2  k  +  l)n 


fc  =  0, 1, JV. 


(31) 


L2(JV  +  1)J’ 

Substituting  Eq.  (28)  into  Eq.  (30)  produces  a  system  of  linear  equations  for  the  unknown 
expansion  coefficients;  namely, 


N  - 

J2 <[A(1  +  Vk)Tj(Vk)  ~  —Mvk)}  =  A(1  +  Vk)F(Vk),  k  =  0, 1, N.  (32) 

3= 0 

This  linear  system  can  solved  by  conventional  means  without  difficulty.  For  sake  of  com- 
parision  with  the  published  results  of  Thynell  and  Ozisik  [7],  we  need  to  introduce  the 
property  of  reflectivity.  Here,  the  reflectivity  is  defined  as 

Ref  =  27t  f  fjLe~2ttRdjj.+ 

J  n= 0 
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[  (l  +  tj)2'$!(ti)  f  cosh[X(l  +  rj)t]e  >2(i+n)2(i  ^dtdrj.  (33) 
7.77=—  1  7  t=o 

To  illustrate  the  effectiveness  of  orthogonal  collocation,  Table  1  presents  a  suite  of  results 
for  the  reflectivity  as  defined  by  Eq.  (33).  The  results  presented  in  Table  1  use  various 
values  of  u  and  radii  (R).  These  simulations  axe  compared  to  the  results  of  Thynell  and 
Ozisik  [7]_.  Thynell  and  Ozisik  [7]  developed  benchmark  results  using  a  Galerkin  method. 
The  last,  column  in  Table  1  indicates  their  results  and  the  number  of  terms  required  to 
obtain  five  digit  accuracy  in  the  reflectivity  calculation.  It  is  clear  from  viewing  this  table, 
that  the  proposed  Chebyshev  collocation  method  produces  benchmark  accuracy.  Also,  a 
collocation  method  is  computationally  more  efficient  than  a  Galerkin  method. 

Table  1.  Numerical  results  for  the  reflectivity  using  the  proposed  Chebyshev  collocation 
method  for  various  ui  and  R. 


CO 

R 

N=4 

N=6 

N=8 

N=10 

N=16 

RefK/l 

0.1 

0.1 

0.88774 

0.88774 

0.88774 

0.5 

0.1 

0.93580 

0.93580 

0.93580 

0.93580  (2)  | 

0.9 

0.1 

0.98677 

0.98677 

0.98677 

0.1 

1 

0.33692 

0.33688^ 

0.33687 

0.33687 

0.33687  (5) 

0.5 

1 

0.54243 

0.54230 

0.54228 

0.54227 

0.54226 

0.54226  (6) 

0.9 

l 

0.87809 

0.87803 

0.87802 

0.87801 

0.87801  (4) 

0.1 

10 

0.02904 

0.02978 

0.02966 

0.02956 

0.02951 

0.02948 

0.02947 

BBCTK3I 

0.5 

10 

0.17085 

0.17207 

0.17128 

0.17084 

0.17062 

0.1705 

0.17044 

WSEsBMSBM 

0.9 

10 

0.54837 

0.54790 

0.54718 

0.54676 

0.54656 

0.54645 

0.54639 

EEEESKEB 

1 

2.4  Orthogonal  Collocation  with  Chebyshev  Basis  (Linear  Anisotropic  Case). 

After  carefully  reviewing  several  papers,  the  derivation  offered  by  Thynell  and  Ozisik  [5] 
appears  preferable  to  other  works.  Following  Thynell  and  Ozisik  [5],  the  general  integral 
form  of  the  transport  equation  in  spherical  coordinates  can  be  expressed  as 


rR  N 

rGn(r)  =  2n  J  ^Q(2')^'0n(7'j  d*  ^  )  omGm(x)A:TTin(r,  x) 


xdx+ 


2 irr  f1  [(-l)ner^  +  e-r#i3 r{R,-^0{r,p4)\Pn{pi)e-^R'r^ xe[OtR],  n  =  0, 1, ...,  N, 

7^=o 

(34) 


where  the  kernel  function  kmn(r, x)  is  given  by 


^mn^)  —  [  Pm( 

Jt=\r-x\  V 


r2  —  x2  —  t2 
2xt 


)Pn{ 


r2-x2+t2  -t 


2rt 


)6 

— dt,  m,n  =  0, 1, TV, 
t 


while 

and 


a(x, r,  fj.)  =  yjx*  -  r2(  1  -  At2), 


(35) 

(36) 


Mo(r,  A4)  =  \/l  -  (r/-R)2(l  -  M2),  (37) 

Assuming  a  transparent  boundary  at  r  =  J?  allows  us  to  write  the  boundary  condition  [7] 

I~(R,-p)  =  f(p),  fj.>0.  (38) 

Here,  the  conventional  form  of  the  phase  function  is  assumed;  namely, 

N 

=  Y^amPm(^)Pm(^),  O.0  =  1,  (39) 

771  =  0 

where  the  {om}^_0  are  known  coefficients. 

For  linear  anisotropic  scattering,  N  =  1,  therefore  we  can  express  Eq.  (36)  explicitly  as  a 
system  of  coupled  weakly-singular  Fredholm  integral  equations  of  the  second  kind.  Thus, 
knowing  Po{p)  =  1,  Pi(a0  =  A4  an^  assuming  p{  —  pf  =  0,  T(p)  =  l,  A  =  0,  and  an 
isotropic  external  source,  f(p)  =  /  [5],  we  find 

rR  w  1 

rGo(r)=2n  [Q(x)k00{r,x)  +  —  amGm(x)kmo(r,x)]xdx+ 

J*= o  47r  m^O 

47rr  /  cosh(rp)t 
Jn= o 

and 


0-a(R,r,n) 


fdn, 


(40) 


rGi(r)  =27 r  /  [Q(x)koi(r,x)  +  —  ^  amGm(x)A;ml(r,x)]xdx- 
Jt=0  771=0 

Anr  f  sinh{rp)e~a^R'r'^ f  pdp,  (41) 

Jn=o 

Next,  we  map  Eqs.  (40,41)  into  a  convenient  form  for  the  introduction  of  the  Chebyshev 
basis  functions.  Let 

r  =  A(1  +77),  re[0,  J2],  ^[-l.l],  (42) 

where  A  =  R/2.  Making  use  of  the  independent  variable  transformation  offered  in  Eq. 
(1.4.4),  the  system  of  linear  Fredholm  equations  displayed  in  Eqs.  (1.4.3a,b)  become 

r1  u  1 

A(1  +  77)4*0(77)  =  2ttA2  /  [QiQKuM  +  j~Y,  «m*m(0*mofa0](l  +  SR+ 

J£=-1  q7r  m=n 
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4nfX(l+r])  j  cosh[X(l  +  r])n]e  a(RMi+v),n) ^  (43) 

“  >m=0 

and 

/•l  u  1 

A(l+.77)$1(77)  =  2?rA2  /  [Q(S)Koi(r),€)  +  —  ^  aTO$m(0-Kmi(»7,£)](l +  0<*£- 

,'^=~1  n  m=0 

47r/A(l  +  7])  f  sinh[X(l +r])n]e~0,(~R’x('1+v^1^ /j,dn,  (44) 

Jn= 0 

where 

$m(v)  =Gm(X(l+r]),  (45) 

Q(rj)=Q(X(l+V),  (46) 

Kmn(v,0  =  femn(A(l  +77),  A(1  +  0).  (47) 

The  coupled  Fredholm  integral  equations  are  in  a  proper  form  for  the  introduction  of  the 
Chebyshev  basis.  Let 


*oM  ■£»!  « 

j—  0  j=0 


♦ita)  =  ='EC?°T>W,  (49) 

j= 0  j=0 

where  Tj{rj)  represents  the  jth  Chebyshev  polynomial  of  the  first  kind.  (Note:  it  is  not 
necessary  that  No  =  N\.)  For  convenience,  we  let  Nq  =  Ni  =  N.  Upon  substituting 
the  approximations  shown  in  Eqs.  (34,35)  into  Eqs.  (43,44),  we  form  the  corresponding 
residual  equations  for  the  system  of  interest;  namely, 

,.\2  /  /*l 

RgW  =  -A(l  +  i,)«u)  +  —(00  J  <(«)ifoo(t;.f)(l  +  i)<H+ 

OJ  /__l*f'K)Xio(»I.fl(l +«)<*{)  +90W,  (50) 

and 

wA2  /  /*i 

(v)  —  ~ A(1  +  v)^{v)  d — 2~  \a°  J  *ZWoi(v,S)(l +  $)<%+ 


0,1  J  ^liOKniv,  €){!■+ €)d£,j  +gi(rj), 


9 


where 


90 (v)  =  9o(v) '+  27rA2 

[  Q(0(!  +  €)Koo(t],  £)d4, 

(52) 

9i(v)  =9i(v)  +  2ttA2 

'€=-1 

(53) 

9o(v)  =  4tt/A(1  +  77)  f 

cos/i[A(l  +  77)/i]e-“^'R’A^1+T?)’^dp, 

(54) 

9iiv)  =  -4tt/A(1  +  7?)  f 

Jfi=  0 

sinh[A(l  +  tj)  dp, 

(55) 

where  it  is  understood  that  R^(rj),  j  =  0, 1  represents  the  residual  functions. 

The  collocation  method  is  defined  through  the  following  orthogonality  statement  [11] 

(R?(V),  -  Vk))i  =  0,  k  =  0, 1, ...,  N,  (56) 

where  the  collocation  points  r)k  can  be  defined  (preliminary)  through  the  open  rule 

^  =  0,1, -,1V.  (57) 


2.5  Orthogonal  Collocation  with  Chebyshev  Basis  (Anisotropic  Case  [12]). 

A  benchmark  numerical  study  has  been  identified  for  the  linear  isotropic  scattering  problem 
in  a  sphere.  Most  studies  consider  the  spherical  shell  problem  over  the  solid  sphere  problem. 
A  careful  study  was  found  [12]  in  the  literature.  The  problem  described  by  Thynell  [12] 
involves  a  uniform  source  with  a  zero  boundary  condition  at  the  outer  radius.  Thynell  [12] 
used  the  integral  form  of  the  transport  equation  and  developed  his  solution  based  on  a 
Galerkin  approach.  The  expansions  for  the  zeroth  and  first  moments  of  the  intensity  was 
based  on  polynomials.  It  is  interesting  to  note  that  these  expansions  were  not  directly  used 
in  developing  plots  of  these  functions  or  even  in  obtaining  the  hemispherical  emissivity.  In 
fact,  Thynell  apparently  obtains  results  by  an  additional  smoothing  process  involving  the 
original  Fredholm  integral  equations. 

For  sake  of  concreteness,  we  begin  by  stating  the  assumptions  used  in  [1]  and  reduce  our 
previously  developed  equations.  Let  /  =  0,  Q{r})  =  1  —  u>,  and  a\  =  +2.7, —2.7  (i.e., 
Thynell’s  g'  =  0.9,  —0.9,  where  a\  =  3 g')  allows  us  for  form  the  equivalent  problem  stated 
in  Ref.  [12]  for  comparative  purposes.  Equations  (43,44)  can  be  shown  to  reduce  to 


A(l+77)$0(?7)  =  27rA2 


i 

-w)Koo{rj,£)  +  —  °m^m(s)-^mo(77,s)](l  +  O^S>  (58) 
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and 


pi  1 

A (1+77)$]. (77)  =  27rA2  /  [(l-U;)-Koi(77,0  +  ^;  '^2  °m$m  (O-^ml  fa,  0]  (l+O^O  fa9) 

*'£=-1  m= 0 


Introducing  the  approximations  for  the  dependent  variables;  namely, 


No 


Mv)  =  J2b>Tj(v) »  <"fo)  =  E6f : t-W. 

j=o  i=o 


and 


^1 


=  X>,t+)  «  *?■(>))  =  Ecf  tiW. 

j— 0  i=o 


(60) 


(61) 


where  7}  (77)  represents  the  jth  Chebyshev  polynomial  of  the  first  kind  allows  us  to  formulate 
the  necessary  residual  equations  from  which  the  collocation  method  can  be  applied.  For 
convenience,  we  let  Nq  =  Ni  =  N.  Upon  substituting  the  approximations  shown  in  Eqs. 
(60,61)  into  Eqs.  (58,59),  we  form  the  corresponding  residual  equations  for  the  system  of 
interest;  namely, 

\  2  />l 

Rq{v)  —  — Ml  +Tl)$o(rl)  H — 2~  (a°  J  (O-ffaofa, 0(1  +  0<^+ 

Oi  f  $f(0#iofa,0(l  +  0ds)  +  Pofa)> 

Ji=-i  J 


(62) 


where 


R? (t?)  =  -A(l  +  r>)*?(ri)  +  («o  f  <(0*oifa,  Oft  +  04c+ 

J  i*iV(0-K’ii(77»0(l+0^)  +51(7), 

50(77)  =  27rA2(l-u;)  f1  (1+0^00(77,0^, 

^=-1 

Pi (77)  =  2ttA2(1  - w)  [  (1  +  O^oi(77, 0^, 

Js=-l 

since  50(77)  =  Pi (77)  =  0  (i.e.,  owing  to  /  =  0). 

The  collocation  method  is  defined  through  the  following  orthogonality  statement  [11] 

(77), £(77  ~  T7fc))1  =  fc  =  0,  l,...,iV,  f  =  0,l, 


(63) 

(64) 

(65) 


(66) 
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where  the  collocation  points  rjk  can  be  defined  through  the  open  rule 


(2k  +  1)tt\ 

2  (N  +  1))' 


k  =  0,1, 


N. 


(67) 


Explicitly  stated  in  terms  of  the  expansions,  the  resulting  system  for  the  unknown  expan¬ 
sion  coefficients  becomes 


£>£[-A(i+,OTmta)  +  ^ f_  (l  +  f)Tm({)Jfooht,«K]+ 

•m  — D  ^  ^ 


E  £jl  +  ()Tm(()K10Md4}=  -sofa*), 


(68) 


and 


f>£[ ^ (i+()Tmmai{rik,m\ 
tZo  1  1  ■/«=-1 


+ 


£i[- A(1  /'_  (1  +«T„({)A-1i(Tft,«)<Je]=  -Si (It), 

m=0  1 


k  =  0,l,...,N. 


(69) 


A  system  of  linear  algebraic  equations  containing  2iV  +  2  unknowns  for  the  expansion 
coefficients  is  obtained  from  the  collocation  method.  The  major  computational  effort  of 
this  method  lies  in  arriving  at  the  numerical  values  for  the  integrals  displayed  in  Eqs. 
(68,69).  Once  the  expansion  coefficients  are  determined,  reconstruction  of  $q  (v)  and 
$£^(77)  is  obtained  through  Eqs.  (60,61),  respectively. 


Tables  2  and  3  present 
defined  by  [12] 


numerical  results  for  the  so-called  hemispherical  emissivity  as 


*i(l) 


(70) 


for  two  optical  thicknesses  (R  =  1,5)  where  A  =  R/ 2,  two  linearly  anisotropic  scattering 
phase  functions  (ax  =  -2.7, 2.7)  and  for  three  distinct  albedos  u.  These  tables  also  indicate 
the  number  of  terms  N  retained  in  the  finite  series  representations  for  (77)  and  (77). 
Here  ox  =  3 g'  where  g'  is  used  by  Thynell  [12].  As  noted  by  Thynell  [12]  and  observed 
here,  convergence  to  an  accurate  spatial  distribution  for  both  (77)  and  (77)  is  rapid 
and  fairly  independent  of  R,u,ai .  It  should  be  noted  that  the  mathematical  description 
for  the  phase  function  can  drastically  change  the  spatial  distribution  of  both  $q  (77)  and 

^(77). 

Tables  4  through  7  present  the  expansion  coefficients,  b f  and  required  for  obtaining 
(77)  and  (77)  when  R  =  1,5;  w  =  0.5,  and  ax  =  -2.7.  It  should  be  noted  that  the  basis 
functions  are  bounded  by  unity,  i.e.,  |Xy  (77)  |  <  1  and  thus  the  significance  of  the  expansion 
coefficients  required  in  developing  (v)  (77)  can  be  evaluated.  The  coupling  effect 


12 


among  the  coefficients  is  evident.  Also,  qualitative  convergence  in  the  numerical  value  for 
the  expansion  coefficients  become  clear  as  N  is  increased. 

Figures  1  and  2  present  the  spatial  distributions  of  &q  (77)  and  $^(77)  when  u  =  0.5, 
ai  =  —2.7  and  R  =  1,5,  respectively.  It  is  evident  that  graphical  accuracy  is  established 
rapidly.  Graphical  convergence  to  accurate  solutions  for  $q  (77)  and  ${^(77)  is  reached  when 
N  —  8  for  R  =  1  and  N  =  10  for  R  =  5.  It  should  be  noted  that  the  majority  of  the 
computational  effort  needed  in  the  proposed  procedure  lies  in  performing  the  integrations 
displayed  in  Eqs.  (68,69). 


Table  2.  Comparison  of  results  for  the  emissivity,  e  using  the  proposed  method  with 
that  of  a  Galerkin  solution  [12]  when  R  —  1  for  various  u  and  N. 


0 

N 

a,  =  -2.7 

O) 

II 

K) 

-4 

0.1 

4 

0.656129 

0.665831 

8 

0.658149 

0.667812 

12 

0.658311 

0.667973 

Exact . 

0.658359 

0.668021 

0.5 

4 

0.444454 

0.469481 

8 

0.445653 

0.470625 

12 

0.445756 

0.470728 

Exact  ^ 

0.445788 

0.470760 

0.9 

4 

0.11997 

0.123487 

8 

0.120238 

0.123775 

12 

0.120264 

0.123781 

Exact . 

0.120273 

0.123789 
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Table  3.  Comparison  of  results  for  the  emissivity,  e  using  the  proposed  method  with 
that  of  a  Galerkin  solution  [12]  when  R  =  5  for  various  u  and  N . 


© 

N 

r- 

II 

nT 

a,  =  2.7 

0.1 

5 

0.907208 

0.947699 

10 

0.93191 

0.970297 

15 

0.933751 

0.972101 

20 

0.934103 

0.97245! 

Exact 

0.934284 

0.972631 

0.5 

5 

0.70546 

0.890032 

10 

0.721914 

0.90223 

15 

0.723144 

0.903432 

20 

0.72339 

0.903681 

Exact 

0.723522 

0.903817 

0.9 

5 

0.323634 

0.455971 

10 

0.327394 

0.460083 

15 

0.327792 

0.460546 

20 

0.327877 

0.460646 

Exact 

0.327928 

0.460706 

Table  4.  Expansion  coefficients  for  QoiVi)  when  R  =  1,  ax  =  —2.7,  and  u>  =  0.5. 


j 

N  =  4 

00 

II 

* 

II 

to 

0 

4.86016 

4.85804 

4.85771 

1 

-1.72849 

-1.73392 

-1.73467 

2 

-0.596946 

-0.603895 

-0.60471 

3 

-0.123983 

-0.133823 

-0.134752 

4 

-0.0403582 

-0.056169 

-0.0572673 

5 

-0.0268265 

-0.0281706 

6 

-0.0144178 

-0.0161221 

7 

-0.00770068 

-0.00992538 

8 

-0.00337749 

-0.0063966 

9 

-0.00419902 

10 

-0.00271823 

ii 

-0.00163689 

12 

-0.000768877 
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Table  5.  Expansion  coefficients  for  $1(77)  when  R  =  l,  ai  =  —2.7,  and  u  =  0.5. 


j 

II 

V 

00 

11 

cf ,  N  -  12 

0 

0.633728  ' 

0.633851 

0.633855 

1 

0.683529 

0.683839 

0.68385 

2 

0.0628458 

0.0632614 

0.0632769 

3 

0.0145621 

0.0152693 

0.015288 

4 

0.00162813 

0.00261631 

0.00263845 

5 

0.000831802 

0.000861553 

6 

0.000294054 

0.000331084 

7 

0.0000954815 

0.000151107 

8 

-0.000000463261 

0.0000746734 

9 

0.0000375295 

10 

0.0000176865 

11 

0.00000463922 

12 

-0.00000441205 

Table  6.  Expansion  coefficients  for  when  R  =  5,  aj.  =  —2.7,  and  u  —  0.5. 


j 

5 

O 

fH 

II 

N  =  20 

0 

10.6347 

10.6287 

10.6274 

1 

-2.94051 

-2.95386 

-2.95645 

2 

-1.73385 

-1.75132 

-1.75414 

3 

-0.831056 

-0.85626 

-0.859423 

4 

-0.364259 

-0.399347 

-0.402927 

5 

-0.131184 

-0.190301 

•0.19442 

6 

-0.0989 

-0.103734 

7 

-0.0557601 

-0.0615263 

8 

-0.0327512 

-0.0398129 

9 

-0.0186327 

-0.0273765 

10 

-0.00836338 

-0.0196322 

11 

-0.0144993 

12 

-0.0109325 

13 

-0.00835725 

14 

-0.00643439 

15 

-0.00495419 

16 

-0.00377851 

17 

-0.00281609 

18 

-0.00199941 

19 

-0.00128297 

20 

-0.00062252 
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Table  7.  Expansion  coefficients  for  $1(77)  when  R  =  5,  ai  =  —2.7,  and  u  =  0.5. 


j 

cf,  N  =  5 

A 

* 

II 

O 

•fW 

* 

II 

8 

0 

0.550292 

0.551624 

0.55166 

1 

0.876375 

0.879222 

0.879314 

2 

0.480069 

0.484099 

0.484231 

3 

0.21294 

0.217996 
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7 

0.00476263 

0.00514869 

8 

0.00198648 

0.00245649 
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0.00130416 

10 

-0.0000276501 

0.000749993 

11 

0.000456405 

12 

0.000289022 

13 

0.000187398 

14 

0.000122979 

15 

0.0000795706 

16 

0.0000496023 

17 

0.0000268658 

18 

0.00000958812 

19 

-0.0000050528 

20 

-0.0000176013 

16 


Figure  1.  §o(v)  and  $1(77)  distributions  for  various  values  of  N  when  u  =  0.5, 
<Zi  =  —2.7  and  R  =  1. 
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3 


Figure  2.  $0 (77)  and  $1(77)  distributions  for  various  values  of  N  when  u  =  0.5, 
ai  =  —2.7  and  R  =  5. 
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Section  3 

Symbolic-Manipulation  of  the  Kernel  Functions 

3.1  Identification  of  the  Kernel  in  Terms  of  Known  Functions. 

In  this  chapter,  we  document  how  to  express  an  highly  anisotropic  scattering  phase  func¬ 
tion  in  terms  of  exponential  integral  functions  and  incomplete  Gamma  functions  by  means 
of  symbolic  manipulation.  This  key  point  now  permits  the  integral  form  of  the  transport 
equation  to  be  used  in  benchmarking  purely  numerical  solutions.  A  major  objective  in¬ 
volves  expressing  the  kernel  function  in  a  way  that  allows  for  the  accurate  integration  of 
the  angular  variable  to  be  performed.  This  kernel  Kmn{r},€ )  was  previously  given  in  terms 
of  a  complicated  integral  involving  products  of  Legendre  polynomials  having  nontrivial 
arguments.  The  ability  to  represent  the  kernel  in  this  new  analytically  integrated  form  is 
of  paramount  importance  for  benchmarking. 

3.2  General  Anisotropic  Scattering  Cases. 

Previously,  we  expressed  the  integral  form  of  the  one-speed  transport  equation  as 


R  N 

rGn(r)  =  J  jc?(x)fcon(r5z)  +  —  ]P  amGm(x)/cmn(r,  x) 


xdx+ 


2tt r  [l  [(-l)ne^  +  e-^]J-[i?,-/i0(r,/z)]-Pn(M)e"a(R’r’M)^, 

J  fj.—0 

xe[0,R},  n  =  0, 1, ...,  N,  (71) 

where  the  kernel  function  kmn(r,x)  is  given  by 


rT+x  /r>2  _  _ +2  \  / rp 2  —  «.2  _i_  +2  v  g  t 

k™(T-X)  =  jL*  P">(  2 Tt  ")P-(  2rt  )~Tdt’  m’n= 


while 

and 


a(x,r,fi)  =  \[xfi  -  r2(l  -  ^2), 


(72) 

(73) 


Po{r,fi)  —  n/T  -  (r/-R)2(l  -  C74) 

Assuming  a  transparent  boundary  at  r  —  R  allows  us  to  write  the  boundary  condition 

I~(R,-p)  =  f(p),  p>  0.  (75) 


Here,  the  conventional  form  of  the  phase  function  is  assumed;  namely, 

N 

p{p,p!)  =  ^P  amPmifTjPmip')-!  a0  =  T  C7®) 

m= 0 


where  the  {am}£=0  are  known  coefficients  and  Pm(p)  represents  the  mth  Legendre  poly¬ 
nomial  of  the  first  kind. 
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Next,  we  map  Eqs.  (71-74)  into  a  convenient  form  for  the  introduction  of  the  Chebyshev 
basis  functions'.  Let  -  - 

r  =  A(1  +rj),  re[0,  R],  1],  (77) 

where  A  =  R/2.  Malcing  use  of  the  independent  variable  transformation  defined  in  Eq. 

(77) ,  the  system  of  linear  Fredholm  equations  displayed  in  Eq.  (71-74)  becomes 
\ 


rl  N 

A(1  +T))$n(ri)  =  27rA2  jQ(O^0«(»7,4)  +  —  arn$m(€)Kmn{V,t)  (1  + 


2ttA(1  +  t i)  l  f(jj)Pn(p)[(- 1)  V(1+7?)/i  +  e-*(l+»»)M]c-«(*,A<l 
Jfi= 0 


n  =  0, 1, IV, 

where 

^miv)  =  Cm(A(l  +  77), 

Q(v)  =  Q(A(1  +  v)i 

Kmn  (t?,  s)  =  kmn (A(l  d*  7]),  A(1  +  £)),  (77,  £)c[— 1,  1], 

with  the  kernel  function  as 


(78) 

(79) 

(80) 
(81) 


KmniVi  0  — 


L 


M2+t?+$) 
A|t7— Cl 


,A2(l+,)2-A2(i  +  02-t2 

/A2(l  +T7)2  -  A2(1  +  02  +  «2\ 

,e  * 

V  2A(1  +£)t  J  n 

l  2A(1  +r})t  ) 

1  t 

m,n  =  0, 1, ...,  N,  (77,  $)€[-!,!]. 


dt, 

(82) 


Equation  (78)  represents  a  system  of  N+ 1  linear  Fredholm  integral  equations  of  the  second 
kind.  These  equations  are  in  a  proper  form  for  the  introduction  of  the  Chebyshev  basis. 

S.S  Identification  of  Kmn(rj,£)  in  Terms  of  Known  Functions. 

The  kernel  function  shown  in  Eq.  (82)  in  terms  of  exponential  integral  functions  and 
incomplete  Gamma  functions  with  the  aid  of  symbolic  manipulation.  That  is,  we  have 
identified  and  expressed  the  product  of  the  Legendre  polynomials  shown  in  Eq.  (82)  as 


m -bn 

£  J=T”(' 7,?)!’  = 

t=— (m-bn) 


/A2(l+,)2-A2(l  +  02-72l 

7A2(l  +  !7)2-A2(l+£)2  +  !2\ 

V  2A(l  +  0<  ' 

1  nv  2A(1  +  rf)t  ) 

m,n  =  0, 1,  ...,1V.  (83) 
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This  rather  complicated  RHS  is  now  expressed  in  terms  of  positive  and  negative  integer 
exponents  in  t  having  a  variable  coefficient  set  denoted  by  F™n  (77,  £).  These  variable 
coefficients  can  be  determined  through  symbolic  computation.  Upon  substituting  Eq.  (83) 
into  Eq.  (82),  we  find 


m+n  a(2+T7+?) 

Kmn(v,0=  E  /  f-'e-'dt.  (84) 

i=-(m+n )  Jt=\\r}-S\ 


This  new  form  is  quite  amenable  to  further  interpretation.  Before  proceding  further,  we 
state  two  integral  relations  that  will  be  used.  It  can  be  shown  that 


r  «-*  *  Ep(z)  E„(y)  „ 

y,„ V* "  1FT " -3=T-  P-1'2’-’ 


yp- 


and 


f  tp  xe  tdt  =  T{p,z)  -T(p,y),  p  =  1,2,... 

Jt-Z 


(85) 


(86) 


Here.  Ev{u)  represents  the  exponential  integral  function  oipth  order  and  T(p,  u)  represents 
the  pth  incomplete  Gamma  function.  For  convenience,  we  define 


z(v,0  = 

y{v,€)  =  A(2-h?7  +  ^), 

thus,  the  kernel  function  shown  in  Eq.  (84)  may  be  written  as 


(87) 

(88) 


0 

ry(v,£) 

=  E  *r*fo.e) 

1  t*  ^e  *  dt-\- 

(m+n) 

't=z(r?4) 

ry(v£) 

£)  /  f~1e~tdt, 

m,n  =  0, 1,..,  IV. 

Jt=z{Tj£) 

m+n 


In  the  first  summation,  replace  —i  by  i  to  get 

m+n  r$(i 7,0  r-t 


(89) 


ryyn&)  e-t 

Kmn(V,Z)  =  E  /  7i+idt+ 


m+n 


"try  /•yfoO 

E^mn^^)/  m,n  =  0,1,..,  IV,  (i7,4M-l,l]-  (90) 

i=l  Jt=z(v,i) 


Using  the  integrals  shown  in  Eqs.  (85,86),  we  observe  that  the  kernel  can  be 

expressed  as 
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m+n 


Kmn(l 1,0=  XT 


t= 0 


fEj+1(z(r},$)) 


Ei+i(v(v,£))\  . 
yHv,0  ' 


m~hn 

1=1 

This  representation  is  correct  but  does  appear  to  cause  some  numerical  difficulties  when 
m,  n  are  greater  than  or  equal  to  2.  It  was  found  that  the  evaluations  using  Mathematica™ 
version  2.2  took  an  excessive  amount  of  CPU  time  on  a  SGI  and  that  some  numerical  in¬ 
accuracies  appeared  in  the  higher  moments.  Scrutinizing  this  occurrence  permitted  us  to 
make  several  observations  including:  (1)  the  geometric  singularity  caused  problems  near 
q  =  —  1  (i.e.,  r  =  0);  and,  (2)  the  strength  of  the  singularity  in  the  kernel  functions  at 
q  =  £  grows  with  the  degree  of  anisotropy.  The  first  identified  problem  can  be  corrected 
by  modifying  (tailoring)  the  trial  functions.  The  second  problem  has  been  partially  ad¬ 
dressed.  However,  the  second  problem  has  been  substantially  reduced  and  we  are  presently 
addressing  the  region  about  rj  =  £  in  the  kernel  function  Kmn{?hQ  a  focused  manner. 

The  highly  successful  results  involving  anisotropic  scatter  in  the  slab  geometry  and  our 
previously  reported  work  involving  linear  anisotropic  scattering  in  spheres  do  not  imme¬ 
diately  translate  to  the  corresponding  spherical  case  when  n>2.  Thus,  a  major  effort  of 
the  research  program  involved  performing  extremely  difficult  integrations  in  an  accurate 
and  efficient  manner.  Once  this  is  resolved,  the  computational  procedure  is  quite  clear. 

3.4  Modified  Kernel  Approach. 

In  order  to  reduce  the  observed  numerical  difficulties  associated  with  our  previous  deriva¬ 
tion,  we  offer  the  following  re-packaging  concept.  We  express  Eq.  (82)  as 

Kmn(V,t)=Kmn(rl 7,  <£)± 

pA^^)p^WTv)Y-Tdt’  m-n“°'1--JVl  (,l0e(-1’11’ 

(92) 

and  regroup  with  the  hopes  of  reducing  the  problem  in  the  vicinity  about  77  =  £.  Regroup¬ 
ing  terms  in  Eq.  (92)  permits  us  to  express  Kmn(r),£)  as 


yA(2+77+£) 

Jt=\\v- Cl 


)]).  m,n  =  0, 1, ...,  N,  (tj,$ )e[-l,l].  (91) 


where 


Kmn  (jl ?  0  K-mri  (j)  j  £)  H”  {Vi  0  ? 

K-mniV't  s) 

^  r  /A^I+t?)2  — A2(l  +  Q2  — 12\  /A2(l  +  77)2-A2(l+Q2+t2 

2A(1  +6)t  )  n\  2A(l+77)t 


(93) 


riiiv 

Jt-z{r)£) 


22 


and 


m'"  =  0'1 . *  (94) 


ry(Vi 0  /  _ £  v  /  £  \  6  ^ 

JL(7?;o  Pm  ( 2A(1  +  77) ) Pn ( 2A(1  +  77) )  ~Tdt ’  m,n  =  0,1,-’iV’  (^Oe[-l,lj.  (95) 

These  expressions  can  be  simplified  in  a  similar  manner  as  previously  discussed;  namely, 
we  find 

m+"  ?%(»())  «+!«(>!, 0)1 


KmniVi  £)  — 
t  '  '~l 


and 


KmnM  =  E  ft’n)(v,o[ 


#*(»?»  6 


+ 


m-f-n 


E  /^m,n) (77, f ) I T(i, 2(77, 0)  -  r(i,y(r7,C)) 


1=1 


tfmnfaf)  =^Om’n)(,7)[^l(z(77,0)  -^l(y(77,0)]  + 

m+n 

E  9im,n)  (v)\T(i,  zfaO)  -  T{i,y(ri,Z))], 


i= 1 


(96) 


(97) 


where  fjm,n\r],v)  =  0,  for  all  m,n.  Again,  these  have  similar  features  to 

(777.  7l} 

those  given  by  our  original  F±  ’  (v,£)  with  the  exception  of  the  collapsing  behaviour 
near  7?  =  £. 


This  new  kernel  representation  behaves  much  better,  however  there  still  exists  some  difficul¬ 
ties.  Asymptotic  expansions  were  attempted  to  alleviate  some  of  the  numerical  difficulties. 
However,  they  did  not  render  enough  accuracy.  In  the  last  month  of  the  program,  a  litera¬ 
ture  survey  revealed  that  some  recent  work  in  the  area  of  Finite  Integral  Part  Integration 
(in  the  sense  of  Hadamard)  may  assist  in  resolving  the  final  dilemma.  Unfortunately,  the 
contractual  period  has  closed  and  the  financial  resources  are  not  available  to  continue  the 
project. 


S.5  Identified  Anisotropic  Scattering  Test  Phase  Function. 

The  traditional  scattering  phase  function  is  expressed  in  terms  of  Legendre  polynomials  of 
the  first  kind;  namely. 

TO 

Pn(/b/x')  =  YlajPj^Pj^'^  (98) 

j=0 

where  {aj}^=Q  represent  known  coefficients.  For  example,  if  00  =  1,  ai  =  —0.56524, 
a2  =  0.29783,  03  =  0.08571,  04  =  0.01003,  05  =  0.00063,  then  the  phase  function  can 
be  graphically  displayed  as  shown  in  Figure  3.  This  particular  phase  function  has  been 
used  often.  It  should  be  noted  that  by  changing  the  last  coefficient  to  05  =  0.800063  a 
drastic  change  occurs  in  the  resulting  phase  function  as  shown  in  Figure  4.  This  test  phase 
function  appears  to  offer  several  features  worthy  in  a  simulation. 
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Figure  3.  Simple  backscattering  phase  function. 


Figure  4.  Highly  lobular  phase  function. 
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Section  4 

Conclusions  and  Recommendations 

The  last  significant  technical  issue  involved  handling,  in  a  general  manner,  the  numerical 
singularity  that  appears  at  the  origin  as  the  degree  of  anisotropy  is  increased.  In  a  general 
sense,  this  is  a  rather  substantial  problem  that  apparently  has  not  been  addressed  in 
the  literature.  For  example,  though  Thynell  and  Ozisik  [5]  developed  the  general  one¬ 
dimensional  integral  form  of  the  transport  equation  for  the  solid  sphere,  they  do  not 
identify  any  potential  numerical  problems  that  could  occur.  In  a  follow-up  paper  by 
Thynell  [13],  the  author  limits  his  discussion  to  the  trivial  problem  of  linear-anisotropic 
scattering.  Likewise,  Abulwafa  [6]  investigates  the  case  of  highly  anisotropic  scattering 
by  recasting  it  into  a  linear- anisotropic  form  (an  approximation).  In  this  way,  the  author 
could  obtain  numerical  results. 

In  the  planar  geometry  (a  nonsingular  geometry),  our  symbolically  enhanced  methodology 
has  been  demonstrated  to  work  extremely  well  (see  [4]).  This  was  again  indicated  during 
the  first  year  of  this  two-year  program.  In  the  case  of  isotropic  and  linear-anisotropic 
scattering,  again  the  approach  has  been  demonstrated  to  work  extremely  well.  This  was 
demonstrated  during  the  first  half  of  the  second  year  of  the  two-year  program.  The  last 
significant  task  involved  the  final  generalization  to  highly  anisotropic  scattering.  A  general 
and  novel  expression  has  been  derived  indicating  that  a  severe  numerical  singularity  is 
present  at  the  origin  for  a  solid  sphere.  A  significant  effort  has  been  underway  to  resolve 
this  final  technical  issue.  As  of  this  report,  this  final  obstacle  has  not  been  overcome. 
This  last  hurdle  appears  resolvable  and  may  have  been  considered  by  other  researchers  in 
different  areas. 

The  development  of  analytic  methods  in  fracture  mechanics  and  the  boundary  integral 
method  may  lead  to  the  proper  handling  of  such  integrals.  These  integrals  are  strongly 
singular  in  the  sense  of  Hadamard  [14-17]  and  involve  the  so-called  finite  integral  part. 
Strongly  singular  kernels  of  this  type  have  appeared,  in  recent  years,  in  the  study  of 
hypersingular  integral  equations  (boundary  element  method)  and  applications  involving 
fracture  mechanics.  The  final  clue  on  the  proper  resolution  of  singularities  associated  with 
the  spherical  problem  may  be  assisted  by  investigating  these  other  areas. 

The  major  recommendation  offered  by  this  investigator  involves  providing  resources  to 
understand  how  the  investigation  of  strongly  singular  kernels  associated  with  fracture  me¬ 
chanics  and  hypersingular  integral  equations  can  be  used  to  resolve  the  final  research  issue 
associated  with  this  study.  Once  this  is  resolved,  the  numerical  procedure  for  determining 
the  solution  of  the  resulting  system  of  Fredholm  integral  equations  is  well  established  and 
no  major  hurdle  is  forseen. 
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Summary 

This  final  report  delivers  encouraging  results  indicating  the  merit  of  the  proposed 
computational  methodology.  The  preliminary  stage  involving  the  past  year  verifies  the 
concept  proposed  in  the  original  proposal  to  DNA.  Follow-on  funding  is  now  sought  to 
resolve  transport  questions  involving  spherical  geometries.  Each  proposed  task  of  the 
SOW  is  addressed  in  a  self-contained  manner  and  is  presented  as  a  chapter  in  this  final 
report. 
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0.0  Report  Organization 


This  final  report  is  organized  in  two  major  sections.  The  first  section  involves  ad¬ 
dressing  Task  I  while  the  second  portion  of  the  report  directly  addresses  Task  II  as  stated 
in  the  Statement  of  Work.  For  sake  of  reference,  the  tasks  as  stated  in  the  SOW  are: 

I)  Development  of  the  appropriate  mathematical  formulation  and  demonstration  of  the 
proposed  methodology.  Several  benchmark  cases  exist  and  will  be  considered  as  they 
pertain  to  slab  geometries.  The  solution  methodology  will  take  advantage  of  symbolic  ma¬ 
nipulation  and  a  spectral  based  expansion  method  for  determining  the  unknown  functions 
of  interest. 

II)  Formulation  of  the  integral  form  of  the  Boltzmann  transport  equation  in  spherical 
coordinates  in  the  presence  of  a  highly  anisotropic  phase  function. 

Each  task  is  addressed  separately  and  is  meant  to  be  self  contained. 


1.0  Task  I:  Verification  of  the  Computational  Methodology 

This  portion  of  the  final  report  addresses  Task  I  as  reported  in  the  SOW  and  described 
in  the  Report  Organization. 

1.1  Introduction 

In  a  series  of  recent  articles  [1-3],  a  symbolically  enhanced  methodology  has  been  de¬ 
veloped  for  solving  the  integral  form  of  the  Boltzmann  transport  equation  in  the  presence 
of  anisotropy.  The  integral  form  of  the  transport  equation  has  been  utilized  in  these  past 
one-dimensional  studies  since  a  reduction  in  dimensionality  takes  place.  The  resulting  sys¬ 
tem  of  weakly-singular  Fredholm  integral  equations  of  the  second  kind  arises  coupling  the 
various  moments  of  the  intensity  (radiative  trains  port)  or  neutron  flux  (neutron  transport). 

In  [1],  an  important  computational  attribute  was  identified  pertaining  to  the  integral 
form  of  the  transport  equation.  It  was  observed  that  symbolic  manipulation  can  be  used  to 
automate  the  identification  of  the  kernel  function  needed  to  arrive  at  the  integral  form  of 
the  transport  equation.  Theoretical  considerations  involving  convergence  and  error  anal¬ 
yses  were  established  in  [2,3]  based  on  a  weighted  residual  methodology  where  Chebyshev 
polynomials  were  used  as  the  basis  functions  in  the  expansions.  The  convergence  rate  was 
established  for  the  isotropic  scattering  case  when  implementing  orthogonal  collocation.  In 
[3],  LaClair  and  Frankel  investigated  linear  anisotropic  scattering  and  developed  rigorous 
error  estimates  based  on  functional  analysis  that  extended  the  error  estimates  developed 
in  [2]. 

This  report  presents  the  actual  implementation  of  symbolic  manipulation  for  auto¬ 
matically  setting  up  the  kernel  functions  needed  in  the  integral  form  of  the  transport 
equation  given  an  arbitrary  scattering  phase  function.  Extension  to  two-dimensional  cases 
now  appears  feasible.  Chebyshev  polynomials  of  the  first  kind  are  chosen  as  the  basis  set 
for  approximating  the  unknown  Legendre  moments.  The  unknown  expansion  coefficients 
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axe  obtained  in  a  weighted  residual  sense  using  orthogonal  collocation.  Comparisons  with 
existing  solutions  qualitatively  indicate  the  accuracy  of  the  approach.  This  form  of  er¬ 
ror  evaluation  is  chosen  since  a  functional  analysis  approach  becomes  unwieldy  and  thus 
impractical  for  cases  having  a  high  degree  of  anisotropy. 


1.2  Formulation 

The  general  one-dimensional,  gray  (one-speed),  slab  formulation  follows  from  previous 
presentations  [1,4]  and  thus  is  not  derived  here.  The  independent  variables  are  mapped 
into  a  coordinate  system  conducive  to  the  spatial  format  of  the  spectral  basis  functions  to 
be  shortly  introduced.  Using  r  =  A(1  +  x),  the  derivation  developed  and  described  in  [1,4] 
can  be  expressed  in  the  equivalent  form 

<5„(x)  =  Fn(x)  -(-  —  am  j  Km<n(X\x  —  x0\)Gm(x0)dx0, 

1  m=0  Jxo=-l 


where 


=  0, 1,  ...,M,  xc[— 1,1], 

(1.1a) 

)  =  Rm,n{ A(l  -  X0))<5m,n(A|x  -  X0|), 

(1.16) 

are  given  as 

\\  f  Pm  ( M) M )  -*!*-*•! 

A)  -  j  e  <*  d/1, 

Jn=0  V 

(1.1c) 

T  \\  _  /  x  —  xa  >  0; 

o))  \(-l)m+n,  x0  —  x  >  0,  --  ■ 

(l.ld) 

while  the  mapped  forcing  function  becomes 


Fn{x)  =  2ir  [  Pn(n)[I+(- l,n)ez*Fsl  +  {-l)nr{l,-n)ez^zll]dfi+ 

J  M=0 

27tA  f  s(x0).Ko,n(A|x  -  x0\)dx0,  xe[— 1,1],  (l.le) 

where  A  =  rd/2.  In  [1],  the  nth  Legendre  moment  of  the  intensity  (or  angular  flux)  is 
defined  as  Gn{r)  for  re[0,  rd]  where  rd  is  theoptical  depth.  The  nth  Legendre  polynomial 
of  the  first  kind  is  denoted  as  P»(/i).  Thus,  Gn(x)  =  Gn( A(1  +  x))  represents  the  mapped 
form  of  the  dependent  variable.  Here,  the  coefficients  are  the  prescribed  phase 

function  coefficients. 

Frankel  [1]  indicates  that  symbolic  manipulation  can  be  exploited  for  automating  the 
evaluation  of  the  integral  function  displayed  in  Eq.  (1.1c)  for  ^m,„(A|x  -  x0|)  since  it  is 
readily  recognizable  in  terms  of  a  finite  set  of  exponential  integral  functions.  In  particular, 
the  partial  kernel  function  (?m,n(A|x  —  x0|)  can  be  expressed  as 
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c(m-n)jE;j+1(A|x-x0|)1  (m,n)  =  0, 1, M,  (1.2) 


Smn~  i 

Qm,n(A|l-I0|)  = 

i=0 

where  p^z)  denotes  the  jth  exponential  integral  function  [5].  Given  the  integer  values  of 
77i,  n  from  0  to  M,  the  coefficients  Cjm‘n)  and  number  of  terms  in  the  sum  Nmn ,  can  be 
established  with  the  aid  of  symbolic  manipulation  without  any  difficulty.  In  fact,  Nm,n 
is  known  a  priori  from  properties  of  polynomial  multiplication  [1].  Figure  1  displays  such 
a  cell  of  code  as  prepared  using  the  symbolic  software  package  Mathematica,™  [6]  and 
implemented  on  a  NeXTstation  having  16  MBytes  of  memory.  Similar  functionality  exists 
on  most  conventional  symbolic  manipulators  such  as  Maple  and  MacSyma.  Finally,  it 
should  be  evident  from  viewing  Eq.  (1.1c)  that  =  c^n,m*  and  7Vm  n  =  lVn  m. 

The  product  of  the  Legendre  polynomials  Pm(/i)  and  Pn{^)  is  formed  and  denoted 
as  “fmn”.  Using  the  intrinsic  function,  CoefficientList,  a  list  is  built  that  contains  the 
coefficients  of  the  resulting  multiplication  of  the  two  polynomials.  The  generated  list 
contains  the  coefficients  in  ascending  order  starting  with  the  coefficient  for  /i°  and  ending 
with  the  coefficient  for  /im+n.  The  total  length  of  the  list  is  m  +  n  + 1.  Use  of  the  intrinsic 
function  Length  counts  the  number  of  elements  in  the  generated  list.  Clearly,  this  should 
be  equal  to  m  +  n  +  1.  This  result  is  stored  in  the  array  named  “arraymaxc”  having 
elements  named  “maxcoeffcjn^n]” .  The  coefficients  c('JTn,n^  are  then  constructed  and  stored 
for  later  used.  The  output  portion  of  the  cell  indicates  the  resulting  operations  to  produce 
c(m,n)  fQr  tjie  Rayleigh  scattering  phase  function,  i.e.,  Af  =  2,  a0  =  1,  aj  =  0,  and  a2  =  0.5. 

In  general,  approximately  half  of  the  coefficients  for  j  =  0, 1, ...,  N(m,n)  —  1  are  zero 

in  value  for  fixed  (m,  n). 

1.3  Solution 

A  weighted  residual  methodology,  utilizing  orthogonal  collocation,  has  been  recently 
demonstrated  to  yield  accurate  results  (2,3).  In  this  section,  an  orthogonal  collocation 
method  is  presented  for  finding  an  approximate  solution  to  Gn(x)  as  indicated  in  Eq. 
(1.1a).  Let  the  unknown  function  Gn{x )  be  formally  represented  by  the  series  expansion 

OO 

Gn(x)  =  ^2  bn,iTi(x),  xe[— 1, 1],  n  =  0, 1, ...,  M,  (1.3a) 

»=o 

where  the  basis  functions  {Tt(*7)}So  3X6  chosen  as  the  Chebyshev  polynomials  of  the  first 
kind  [2,3,7]  and  are  expressible  as 

T«(f?)  =  cos[i(cos_177)],  i  =  0, 1, ...  (1.36) 

Chebyshev  polynomials  [7]  have  numerous  exploitable  features  and  have  successfully 
been  used  in  studies  involving  fluid  mechanics  [8],  solid  mechanics  [9],  and  radiative  trans¬ 
port  [10,11].  The  unknown  expansion  coefficients  for  each  moment,  n  requiring  resolution 
are  denoted  as  {6n,»}m=o-  practice,  this  series  representation  is  truncated  after  a  finite 
number  of  terms,  say  at  order  Nn,  which  depends  on  the  particular  moment,  n.  For  sake  of 
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simplicity,  let  TVn  =  TV,  that  is,  each  moment  is  approximated  by  a  polynomial  of  the  same 
order.  Thus,  the  TVth-order  approximation  of  G„(x)  is  denoted  as  G%(x)  and  expressed 

by  '  *>  - 

N 

G„(x)*<S«(x)  =  £</r,(x),  n  =  0, 1 . M,  (1.4) 

»=0 

where  b*t  represents  an  approximation  to  6n,„  respectively  for  each  fixed  n,i  in  the  finite 
set. 

The  series  representation  shown  in  Eq.  (1.4)  for  G* (x),  n  =  0, 1  is  substituted 
into  Eq.  (1.1a)  to  produce 


N 


Rn  (X)  +  bn,iTi(X)  =  £»(*)  +  £  6£.<  £ 

3-0 


t=0 


M 


N 


ttnn-l 


m=0 


i=0 


n  =  0, 1, ...,  M;  xe[— 1, 1],  (1.5) 

where  R*  (x)  represents  the  residual  function  required  to  maintain  the  equality  displayed 
in  Eq.  (1.5a).  The  function  /j+ij(x)  represents  the  result  of  the  analytic  integration 
of  the  product  of  the  ith  Chebyshev  polynomial  of  the  first  kind  and  the  kernel  function 
tfm,n(A|x  -  x0|)  defined  in  Eq.  (1.1b).  That  is,  upon  redefining  Qm,n(A|x  -  x0|)  shown  in 
Eq.  (1.1c)  as  Eq.  (1.2)  and  then  forming  the  definition  of  the  kernel,  tfm,n(A|x  -  x0|),  the 
function  (x)  is  identified  as 

[  -Rm,n(A|x  —  x0j)£'J  +  1(A|x  —  X0|)Tj(x0)cfx0, 

j  =  0,l,...,TVm.„  -  1;  (m,n)  =  0,1,...,  Af;  i  =  0, 1, ..., TV;  xe(-l,l],  (1.6a) 

Integrating  Eq.  (1.6a)  analytically  yields 


#£?(*>  =  E  (-D*£j+*«(0)[i  +  (-!)”>«+*]+ 
0 


E  ’(-UiWWi  + 1))  -  1  -  x))  j , 

j  =  0,l,...,TVm>n  -  1;  (m,n)  =  0,1,..., TW,  i  =  0, 1, ...,TV;  xe[-l,l].  (1.66) 

Here,  the  kth  derivative  of  the  pth  Chebyshev  polynomial  of  the  first  kind  is  denoted  as 

Unless  the  exact  solution  to  Gn(x),  n  =  0,1,...,  is  a  linear  combination  of  the  ba¬ 
sis  functions  {Tt(x)}^0,  no  set  of  coefficients  6^  can  make  R%(x)  vanish  for  all  n  = 
0, 1,...,M,  xe[— 1, 1].  However,  suitable  expansion  coefficients  can  be  obtained  by  making 
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the  residual  functions  indicated  in  Eq.(1.5)  small  in  some  sense.  Let  the  inner  product  of 
two  reed- valued  functions  gi{t)  and  gi[t)  be  defined  as 


{9uh)Wk  =  /  Wi(s)5l(5)p2(s)d5, 

and  the  corresponding  norm  as 


— 1 


lift  II 


■[/ 


wk{s)gl{s)ds 


s=-l 


(1.7  a) 


(1.76) 


where  wk(s)  is  a  non- negative,  real  and  integrable  weight  function. 
For  the  collocation  method,  the  orthogonality  relation  becomes 


(R%{x),5{x  -*k))l  =  0,  n  =  0, 1, ...,  M,  k  =  0, 1, N.  (1.8a) 

Here,  the  Dirac  delta  function  is  denoted  by  6  while  the  N  +  1  collocation  points  are 
indicated  as  xk,  k  =  0, 1, ...,  N  and  can  be  defined  by  the  closed  rule  [12] 


Trfc 

xk  =  cos(— ),  =  0, 1, ...,  N. 


(1.86) 


By  choosing  this  set  of  N  +  1  collocation  points,  the  residual  functions  at  the  endpoints 
of  the  physical  domain  become  zero.  This  is  important  since  key  physical  properties  are 
obtained  from  the  endpoints.  Note  that  from  Eq.  (1.86),  one  interprets  that  i0  —  1  and 
xn  —  —1.  Error  and  convergence  analyses  have  been  performed  illustrating  the  merit  of 
this  choice  of  basis  functions  and  collocation  points  in  the  study  of  the  radiative  equation 
of  transfer  in  both  an  isotropically  [2]  and  a  linear  anisotropically  scattering  medium  [3]. 

Applying  the  orthogonality  condition  displayed  by  Eq.  (1.8a)  to  Eq.  (1.5)  for  k  = 
0,1  produces 


£  <*!■«(**)  C,i  £  e«— >j£»>(„)  .  #„(**), 

i=0  m=0  i=0  j=0 


k  =  0, 1, ...,  iV,  n  =  0, 1,...,M.  (1.9) 

This  linear  system  of  (M+1)(N+1)  equations  for  the  unknown  expansion  coefficients  6^, 
n  =  0, 1, ...,  M,  t  =  0, 1, ...,  N  can  be  solved  by  matrix  inversion.  Reconstruction  of  G%(x) 
can  be  accomplished  with  the  aid  of  the  expansion  shown  in  Eq.  (1.4).  Key  physical 
properties  can  be  deduced  from  the  accurate  knowledge  of  the  zeroth  and  first  moments. 

1.4  Results  and  Discussion 

Numerical  results  are  now  presented  to  illustrate  the  accuracy  of  the  symbolically  en¬ 
hanced  methodology.  It  is  evident  that  the  slab  situation  described  here  is  of  limited  utility. 
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However,  extension  of  such  a  methodology  to  spherical  problems  now  appears  feasible.  All 
of  the  numerical  results  presented  in  this  report  were  developed  using  Mathematica™ 
as  implemented  on  a  NeXTstation.  The  entire  computer  code  developed  for  the  general 
anisotropic  scattering  case  is  less  than  150  lines  in  length.  This  includes  the  graphical  and 
tabular  output  statements.  Theoretical  considerations  involving  convergence  and  rigorous 
error  estimates  have  been  performed  by  Frankel  [2]  and  LaClair  and  Frankel  [3]  for  the 
case  of  isotropic  and  linear-anisotropic  scattering,  respectively. 

Several  “benchmark”  problems  appear  in  the  thermal  radiation  literature.  It  is  evi¬ 
dent  that  the  present  methodology  is  applicable  to  the  neutron  transport  equation.  The 
radiation  literature  is  full  of  benchmark  results  and  thus  proof-of-principle  is  indicated  with 
several  well-established  radiation  problems.  For  this  report,  the  scope  of  the  presentation 
is  limited  to  the  problem  involving  uniformly  incident  radiation  of  unit  intensity  at  x  —  —  1 
(t  =  0)  and  where  no  radiation  is  incident  at  x  =  1  (r  =  rd)  [1,4,13,14].  This  gives  rise 
to  I+{— =  1,  and  J”(l,—  n)  =  0,  n  >  0.  Furthermore,  let  the  internal  source  term, 
s(x0)  be  negligible.  Therefore,  Fn(x)  reduces  to 

F„(x)  =  2ir  [  Pn(n)e~L^±~L dfi,  xe[— 1, 1],  n  =  0, 1, ...,  Af,  (1.10a) 

Jfi=0 

or  upon  making  a  similar  observation  as  noted  in  developing  the  kernel  functions, 


Fn(x)  =  2tt 

i= o 


J°’n) 

S 


£j+2(A(1  +  x)),  xe[— 1,1],  n  =  0, 1, ....  M. 


(1.106) 


Next,  consider  the  two  scattering  phase  functions  shown  in  Figure  2. where  the  phase 
coefficients  are  presented  in  Table  1.  Figure  2a  displays  a  weakly,  back-scattering  phase 
function  while  Fig.  2b  portrays  a  weakly,  forward-scattering  phase  function.  Both  of  these 
cases  are  well  documented  and  thus  serve  as  the  benchmark  cases  for  consideration. 


Using  PF1  (backward-scattering  case),  Table  2  presents  numerical  results  for  the  re¬ 
flectivity  and  transmissivity,  as  defined  by  [1,4] 


R=1.$±nxl.m=ii. 


(l.llo) 


and 


G,(l)  ..Cftl) 

7T  ~  7T  ’ 


(1.116) 


respectively,  for  various  optical  depths,  rd  and  single-scattering  albedos,  u.  This  table  also 
contains  the  results  of  three  previous  investigations.  The  results  of  Frankel  [1],  Thynell 
and  Ozisik  [4],  and  Cengel  and  Ozisik  [13]  are  presented.  Frankel  [1]  used  singularity 
subtraction  and  trapezoidal  integration  to  obtain  numerical  solutions.  Thynell  and  Ozisik 
[4]  developed  an  intricate  solution  based  on  eigenfunction  expansions.  This  methodology 
required  a  separate  analysis  for  the  conservative  case  where  u>  =  1.  No  special  treatment 
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is  required  by  the  present  approach  when  u  =  1.  Finally,  Cengel  and  Ozisik  [13]  developed 
a  Galerkin  solution  in  which  Legendre  polynomials  of  the  first  kind  were  used  in  the  series 
representations.  It  is  evident  that  the  present  scheme  rapidly  converges  to  the  “numerically 
exact*  solution  offered  by  Cengel  and  Ozisik.  Upon  closer  inspection,  the  present  approach 
appears  to  be  converging  faster  than  the  solution  developed  by  Thynell  and  Ozisik  [13], 

Using  PF2  (forward-scattering  case),  Table  3  compares  the  results  of  Frankel  [1], 
and  Menguc  and  Viskanta  [14]  to  the  present  approach.  It  is  evident  from  viewing  this 
table  that  the  flux  results,  defined  as  Q{ x)  =  G?{ x),  produced  by  the  present  method  is 
comparable  to  the  F-N  solution  when  rd  =  1  and  u  =  0.8. 

Symbolic  manipulators  are  commonly  considered  to  be  RAM  intensive.  Figure  3  illus¬ 
trates  the  memory  requirements  for  a  complete  simulation  as  a  function  of  the  approxima¬ 
tion  order  N  for  PF1  when  M  =  5,  u>  =  0.9,  and  rd  =  1.  Memory  demands  remain  minimal 
and  appear  to  grow  linearly  with  increasing  values  of  N.  A  plot  of  the  RAM  requirements 
as  a  function  of  the  degree  of  anisotropy,  Af,  for  fixed  N  reveals  a  similar  trend.  Actual 
run  times  for  generating  numerical  solutions  using  a  symbolic  manipulator  can  become 
excessive  as  M  and  N  grow.  The  majority  of  the  CPU  time  used  in  a  simulation  arises 
from  constructing  the  coefficient  matrix  and  known  vector  indicated  in  Eq.  (1.9).  This 
situation  can  be  remedied  by  developing  a  hybrid  method  that  uses  a  high-level  language 
for  the  actual  evaluations  of  Ti{xk),  I^[nJ(xk)  and  Fn{xk).  CPU  times  then  become  mini, 
mal  owing  to  the  nature  of  the  programming  language.  Combining  symbolic  manipulation 
with  a  high-level  language  would  make  such  problems  routine.  Contemporary  symbolic 
packages,  such  as  Mathematica™ ,  allow  for  this  type  of  interaction. 

In  concluding  this  portion  of  the  report,  it  appears  that  the  symbolically  augmented 
methodology  used  in  concert  with  orthogonal  collocation  having  a  Chebyshev  basis  set 
performs  well  for  anisotropically  scattering  phenomena.  Interfacing  the  present  approach 
with  a  high-level  language  reduces  the  actual  computational  times. 
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1.6  List  of  Tables 


Table  1:  Phase  function  coefficients  denoted  as  PFl  (backward)  and  PF2  (forward)  corre¬ 
sponding  to  Figure  2. 


Table  2:  Comparison  of  results  for  the  reflectivity,  R  and  transmissivity,  T  among  recent 
studies  using  PFl  for  various  single-scattering  albedos,  ui  and  optical  depths,  rd. 

Table  3:  Comparison  of  results  for  the  surface  heat  fluxes  among  recent  studies  using  PF2 
when  £j  =  0.8  and  rd  =  1. 


A-ll 


Term,  m 


0 

1 

2 

3 

4 

5 

6 


Phase  Function  1  (PF1, 
M=5)  Back  Scattering 
Coefficients 

I 

-0.56524 

0.29783 

0.08571 

0.01003 

0.00063 


Phase  Function  2  (PF2, 
M=6)  Forward  Scattering 
Coefficients 

1 

0.643833 

0.554231 

0.103545 

0.010498 

0.000563 

0.000019 


Table  1:  Phase  function  coefficients  denoted  as  PFl  (backward)  and  PF2  (forward)  corre¬ 
sponding  to  Figure  2. 
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N 

Td  = 

0.1 

u  - 

1 

Td  = 

5 

R 

T 

R 

T 

R 

T 

4 

0  0083562 

0  838491 

0  0250434 

0  229076 

00287262 

0  0019101 

6 

0  0083560 

0  838491 

0  0249877 

0  229064 

0  0266237 

0  0019845 

8 

0.0249822 

0  229063 

0  0263159 

0  0019696 

10 

0  0262535 

0  0019656 

12 

0.0262365 

0  0019645 

100  Panels  [1] 

0.0083562 

0.838491 

0.0249929 

0.229064 

0.0264040 

0.0019642 

N“15  [4] 

0.0083545 

0  83849 

0024982 

0.22906 

0.026227 

0.0019640 

“Exact”  [13] 

0.008356 

0.838491 

0.024981 

0.229063 

0.026226 

0.001964 

4 

0.084895 

0895490 

0388626 

0.439743 

0.513838 

0.0486162 

6 

0.0848944 

0.895490 

0.388459 

0.439786 

0.507288 

0.0417609 

8 

0.388442 

0.439791 

0.506409 

0.041738 

10 

0.506237 

0.0417341 

12 

0.506189 

0.0417334 

100  Panels  [1] 

0.0848961 

0.895490 

0  388541 

0.439800 

0.508286 

0.0417555 

N=15  [4] 

0.084873 

0.89552 

0.38847 

0.43970 

0.50616 

0.041720 

“Exact”  [13] 

0  084895 

0.895489 

0.388437 

0.439792 

0.506157 

0.041733 

4 

0.0959217 

0.904078 

0  486903 

0.513097 

0.821002 

0.178998 

6 

0.0959213 

0.904079 

0.486801 

0.513199 

0.819225 

0.180775 

8 

0.486791 

0.513209 

0.819001 

0.180999 

10 

0.818955 

0.181045 

12 

0.818943 

0.181057 

100  Panels  [1] 

0.0959232 

0.904078 

0  486896 

0.513217 

0.821180 

0.181071 

N=15  [4] 

0.095890 

0.90411 

0.48688 

0.51312 

0.81898 

0.18102 

“Exact”  [13] 

0.095922 

0904078 

0  486788 

0513212 

0818934 

0.181066 

Table  2:  Comparison  of  results  for  the  reflectivity,  R  and  transmissivity,  T  among  recent 
studies  using  PF1  for  various  single-scattering  albedos,  u  and  optical  depths,  rd. 
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Method 

CK-1) 

QO) 

Present  Approach: 

N=4 

0.760344 

0.455868 

N=6 

0.760548 

0.455877 

N=8 

0.760568 

0.455879 

N=10 

0.760573 

0.455879 

Frankel  [1] 

Panels3  40 

0.7601053 

0.4559390 

Panels3  80 

0.7604430 

0.4558964 

Panels3 100 

0.7604881 

0.4558906 

F-N  Method  [14] 

N=3 

0.76061 

0.45596 

N=5 

0.76058 

0.45588 

N=9 

0.76057 

0.45588 

Table  3:  Comparison  of  results  for  the  surface  heat  fluxes  among  recent  studies  using  PF2 
when  u)  =  0.8  and  rd  =  1. 
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1.7  List  of  Figures 


Figure  1:  Mathematica™  input/output  cell  displaying  the  construction  of  the  coefficients 

(m.n) 

\ 

Figure  2:  Anisotropic  scattering  phase  functions  denoted  as  (a)  PFl-weakly,  back-scattering 
(M  =  5);  and,  (b)  PF2-weakly,  forward-scattering  (M  =  6). 

Figure  3:  Memory  requirements  as  a  function  of  N  for  fixed  M  —  5. 
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(•Symbolic  determination  of  coefficients  needed  for 
Kq.  (2)  when  the  degree  of  anisotropy,  M  is  specified. 
The  variable  "alist"  contains  each  of  the  desired 
coefficients  in  ascending  order,  j»0,l,..., 
m&xcoef f c  [m,  n]  -1,  for  a  given  <m,n)  .  Example t 
Rayleigh  scatter,  M-2  *) 

Clear [u,m,n,  fmn,coeffc,maxcoeffc,M] 

,  M-2; 

arrayc-Xrray  [coeffc,  {  (M+l)  A2,M+1,M+1),  0] ; 
arraymaxc-Array  [maxcoef  fc,  {M+1,M+1},0]  ; 

Do [Dot 

Clear [alist] ; 

fmn-Expand  [LegendreP  [m,  u]  *LegendreP[n,u]  ,u]  ; 
alist-Coeff  icientList  [£mn,u]  j 

Print ["alist  [m  •  ",m, ",n  »  ",n, ■  ]  ■  ", alist]; 

m&xcoef  fc  [m,n]  "Length  [alist] ; 

Do  [coeffc  [j  ,m,  n]  "Part  [alist,  j  4-1]  , 

{j,0,maxcoeffc [m,n] -1)  J ,  <m,  0,M>] ,  {n,  0,M)] 


alist 

[m 

= 

0,n 

= 

0 

] 

= 

U) 

alist 

[m 

= 

l,n 

= 

0 

1 

= 

(0,  1} 

1 

3 

alist 

[m 

= 

2,n 

S 

0 

1 

= 

{-(-),  o. 

-} 

2 

2 

alist 

[m 

= 

0,  n 

= 

1 

1 

S 

(0,  1} 

alist 

[m 

= 

l,n 

= 

1 

1 

= 

(0,  0.  1} 

1 

3 

alist 

[tn 

= 

2 ,  n 

= 

1 

1 

= 

{0, 

0, 

-} 

2 

2 

1 

3 

alist 

[m 

= 

0,n 

= 

2 

1 

= 

{-(-),  0, 

-} 

2 

2 

1 

3 

alist 

[m 

= 

l,n 

s 

2 

1 

= 

(0, 

0, 

-> 

2 

2 

1 

3 

alist 

[m 

= 

2  ,n 

= 

2 

1 

= 

0,  -( 

-)  , 

0, 

4  2  4 


Figure  1:  Mathematical™  input/output  cell  displaying  the  construction  of  the  coefficients 
c<m’n). 
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Scattering  Function  PF1 


Figure  2:  Anisotropic  scattering  phase  functions  denoted  as  (a)  PFl-weakly,  back-scattering 
(Af  =  5);  and,  (b)  PF2-weakly,  forward-scattering  (M  =  6).  C 
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Approxim 


e  3:  Memory  requirements  as  a 
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2.0  Task  II:  Spherical  Coordinate  Formulation 

This  portion  of  the  final  report  addresses  Task  II  as  reported  in  the  SOW  and  described 
in  the  Report  Organization. 

2.1  Introduction 

The  development  of  the  integral  form  of  the  Boltzmann  transport  equation  in  the 
spherical  coordinate  system  is  sought  since  the  proposed  computational  methodology  takes 
advantage  of  a  reduction  in  dimensionality  associated  with  the  integral  formalism.  It  is 
interesting  to  note  that  researchers  from  both  the  neutron  and  radiative  transport  are¬ 
nas  have  considered  the  integral  formalism.  However,  in  general,  the  communities  have 
consistently  addressed  the  development  of  numerical  methods  for  solving  the  conventional 
integro-differential  form  of  the  transport  equation. 

A  brief  overview  of  the  pertinent  literature  is  now  offered  from  which  an  appropriate 
integral  formalism  is  presented  for  highly  anisotropic  scatter  in  a  one-speed  situation.  The 
radiative  [1-11]  and  neutron  [12-15]  communities  have  considered  spherical  shell  geometries 
in  numerous  past  studies.  Most  situations  have  considered  highly  idealized  scattering 
situations  mainly  involving  isotropic  scattering.  Thynell  and  Ozisik  [8]  also  considered, 
in  the  context  of  radiative  transport,  the  situation  where  the  single-scattering  albedo  was 
spatially  dependent. 

Thynell  and  Ozisik  [9],  Abulwafa  [10],  and  Pomraning  and  Siewart  [11]  derived  the 
integral  form  of  the  transport  equation  in  spherical  form.  Thynell  and  Ozisik  [9]  developed 
an  integral  formalism  for  a  solid  sphere  while  Abulwafa  [10]  developed  an  integral  formalism 
based  on  a  spherical  shell.  It  is  interesting  to  note  that  Thynell  and  Ozisik  [9]  merely 
derived  the  correct  system  of  weakly-singular  Fredholm  integral  equations  of  the  second 
kind.  That  is,  they  did  not  present  any  numerical  results.  Abulwafa  [10]  presented  some 
numerical  results  as  developed  using  a  Galerkin  method  based  on  a  monomial  set  of  basis 
functions.  The  derivation  offered  by  Abulwafa  [10]  used  a  change  of  variables  [15]  to  arrive 
at  the  desired  integral  formulation.  Owing  to  the  general  nature  of  the  forms  presented  by 
Abulwafa  [10]  and  Thynell  and  Ozisik  [9],  it  is  these  forms  that  appear  most  satisfactory 
for  the  present  research  program. 

2.2  Formulation  and  Remarks 

The  conventional  integro-differential  form  of  the  Boltzmann  transport  equation  in  the 
spherical  coordinate  system  is  [9,10,16] 


-"-^(r,^)+^(r,/x)  =  Q(r)+|  amPm{v)  J  _  ^ 

re[0,  J?],  m«[-  1,1],  (2.1a) 

where  i/>(r,  n)  is  the  angular  flux,  a /  is  the  the  mean  number  of  secondaries  per  collision, 
Q(r)  represents  an  isotropic  source  density  function,  r  is  the  distance  variable  measured 
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in  units  of  mean  free  path,  n  is  the  cosine  of  the  angle  between  the  propagating  flux  and 
the  radial  coordinate,  and  R  is  the  radius  of  the  spherical  domain  of  interest.  It  is  not 
difficult  to  generalize  Eq.  (2.1a)  further  to  include  a  spatially  dependent  u>,  i.e.,  u  -4  u{r) 
[9].  The  phase  function  is  normally  expressed  as 

M 

P( A*»  A*  )  =  ^  ]  a7nPm(P-)Pm(^t  )>  &Q  =  1,  (2-16) 

m=0 

where  the  coefficients  {am}£f_0  are  known  and  Pm(n)  represents  the  mth  Legendre  polyno¬ 
mial  of  the  first  kind.  Often  incoming  distributions  are  prescribed  as  boundary  conditions 
at  the  external  surfaces  (15,16]  in  spherical  shells.  For  the  solid  sphere,  we  consider  the 
isotropic  external  condition  at  r  =  R,  namely 


rp(R,n)  =  F,  fx  =  [-1,1].  (2.1c) 

Further  generalizations  on  the  boundary  condition  are  available. 

Following  Sahni  [15]  and  Abulwafa  [10],  we  introduce  the  neutron  path  coordinates 

£  =  rV(l-/<a),  (2.2a) 


V  =  V-r- 

Recalling  the  chain  rule  from  differential  calculus,  we  find 

d_  =  did_  dvd_ 

dr  dr  d£  dr  dr]  ’ 


and 


d_  _  d£8_  dvd_ 

dfj.  dfi  d£  3/i  dr] 


(2.2  6) 


(2.3  a) 


(2.35) 


Introducing  Eq.  (2.3)  into  Eq.  (2.1a)  produces 
d *  " 


m=0 


+  Q{  +  rp),  (2.4a) 


where 

.  v)  =  i’iV?  —====),  (2.46) 

vr  +  T 

and 

Gm{r)=  f  Pm{n)ip{r,n)dn,  m  =  0, 1, ...,  N,  (2.4c) 

J  pz r— 1 

where  we  use  Abulwafa’s  definition  [10]  of  Gm(r).  (Note  that  Thynell  and  Ozisik  [9]  and 
Abulwafa  [10]  have  slightly  different  definitions  of  the  Legendre  moments  of  ip(r,  p)  in  that 
the  factor  of  2ir  is  omitted  in  Abulwafa’s  definition.)  It  is  evident  that  r  =  y/%2  +  rf2  and 
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H  =  y/y/t2  +  r?2  and  that  re[0,  R]  and  ^te[— 1, 1]  are  transformed  into  \/£2  +  r?2e[0,  R]  and 

^>0. 

At  this  point,  two  approaches  for  arriving  at  the  integral  representation  can  be  fol¬ 
lowed.  Thynell  and  Ozisik  [9]  considered  the  forward  and  backward  components  separately 
as  was  done  by  Frankel  [17]  when  considering  the  slab  geometry.  Abulwafa  [10]  and  Sahni 
[15]  followed  a  direct  route  in  arriving  at  the  integral  form.  Using  either  approach,  an 
integrating  factor  is  introduced  and  upon  a  lengthy  but  relatively  straightforward  set  of 
manipulations  [4],  we  arrive  at  a  system  of  weakly-singular  Fredholm  integral  equations  of 
the  second  kind  [10]  for  the  various  Legendre  moments  of  ^(r,  ji),  namely 


rR  N 

rGn(r)  =  -FbKnl(r, b)+  y[Q(y)Kn0{r , y)  +  -  £  amKnm(r, y)Gm(y)]dy , 

Jy=0  i  m-n 


'V=u  “  m=0 

n  =  0,1,...,JV,  re[0,R], 

where  the  kernel  Knm{r,y)  is  given  as 

>-(r+y)  „.2_Ji_*2, 

-vl 


2  „  0t 


(2.5a) 


(2.55) 


The  weakly-singular  (integrable)  nature  of  these  equations  are  exhibited  in  the  kernel  func¬ 
tions.  Under  certain  conditions,  the  first  exponential  integral  function  Ei(s)  is  produced 
which  contains  a  known  singular  point  at  $  =  0.  Technically,  a  logarithmic  singularity 
exists  at  this  point.  The  exponential  integral  functions  Ek{s),  k  =  2,3,...  are  nonsingular 
and  thus  do  not  represent  any  real  problems. 

For  linear  anisotropic  scattering  ( N  =  1),  Eq.  (2.5a)  reduces  to 


rG0(r)  =  —  FbKoi(r,  6)+  (2.6a) 

J  y[<2(y)tfoo(riy)  +  2  (^oo(riy)Go(y)  +  ai^oi(riy)Gi(y)jjdy, 


and 


rG\{r)  = —FbKn(r,b)+  (2.66) 

JR 

J  y^Q{y)Kio(r,y)  +  ^(K10(r,y)Go(y)  +  a1Kn{r,y)Gi{y)^dy,  re[0,Jl]. 

For  the  linear  anisotropic  case  [10],  the  kernel  functions  can  be  identified  in  terms  of 
exponential  integral  functions.  That  is,  the  kernels  expressed  in  Eq.  (2.6a, b)  are 


-K'oo(riy)  =  Ex (z)  -  Ei(x), 

K0i{r,y )  =  [ssn(r  -  y)yE2{z)  -  Ez(z)  +  yE2{x)  +  Ez{x)]/y, 


(2-7  a) 
(2.76) 
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and 


Kio(r,y)  =  -K01(y,r), 


(2-7  c) 


Ku{r,  y)  =  [Ei(x)  +  ryEz(x)  +  xEA(x)  -  E5(z)  -  zE4(z)  +  ryE3(z)]/(ty),  (2.7 d) 

where  z  =  |r  —  y|  and  i  —  r  +  y. 

\ 

With  this  observation,  the  approach  developed  in  Part  1  of  this  report  appears  possi¬ 
ble.  In  Phase  II  of  the  offered  program,  it  is  our  intent  to  generalize  the  concept  put  forth 
in  Part  1  of  this  report  to  the  spherical  problem  derived  in  Part  2.  Initial  consideration 
will  involve  using  the  proposed  approach  and  confirming  the  results  of  Abulwafa  [10]  who 
has  presented  benchmark  data  for  the  spherical  problem.  Extension  to  highly  anisotropic 
scattering  will  then  be  considered. 
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